function [deltaHedgeEuler, deltaHedgeEulerAV, deltaHedgeEulerD, ...
   deltaHedgeIntegral, deltaHedgeIntegralAV, deltaHedgeIntegralD] ...
   = Digital(steps)

interestRate    = 0.06;
priceZero       = 100;
sigma           = 0.2;
T               = 1;
deltat          = 1 / 200;
N               = T / deltat;
strike          = 99;
epsilon         = 0.9;

iter = 0;
for M = steps
    % Keep track at which position in the result vector the current app-
    % roximation should be stored. This is necessary since the number of
    % steps do not need to be continuous.
   iter = iter + 1;
    % Initalize the Random Number Generator (RNG).
   randn('state', 1);

   for j = 1 : M
       % Generate a random walk phi_t ~ N(0, 1), t <= T
       phi     = randn(N);
       % Generate the random variable for the anti variate.
       phiAV   = -phi;

       % Generate a stock price.
       inter     = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero);
       interAV	  = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phiAV, priceZero);
       interDH      = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero + epsilon);
       interDHAV    = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phiAV, priceZero + epsilon);
       interDHS      = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero - epsilon);


       euler     = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero);
       eulerAV   = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phiAV, priceZero);
       eulerDH   = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero + epsilon);
       eulerDHAV = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phiAV, priceZero + epsilon);
       eulerDHS  = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero - epsilon);

       % Compute the payoff of the two call options. The long has a pay off
       % equal to (S - (K - 0.5))^+ and the short -(S - (K + 0.5))^+ The sum of
       % the two form the portfolio which can be used to replicate an
       % approximation of a digital option.
       valueIntShort    = inter(N + 1) - strike - 0.5;
       valueIntAVShort  = interAV(N + 1) - strike - 0.5;
       valueIntLong     = inter(N + 1) - strike + 0.5;
       valueIntAVLong   = interAV(N + 1) - strike + 0.5;
       valueIntDHShort    = interDH(N + 1) - strike - 0.5;
       valueIntDHAVShort  = interDHAV(N + 1) - strike - 0.5;
       valueIntDHSShort   = interDHS(N + 1) - strike - 0.5;
       valueIntDHLong     = interDH(N + 1) - strike + 0.5;
       valueIntDHAVLong   = interDHAV(N + 1) - strike + 0.5;
       valueIntDHSLong    = interDHS(N + 1) - strike + 0.5;

       valueEulShort    = euler(N + 1) - strike - 0.5;
       valueEulAVShort  = eulerAV(N + 1) - strike - 0.5;
       valueEulLong     = euler(N + 1) - strike + 0.5;
       valueEulAVLong   = eulerAV(N + 1) - strike + 0.5;
	   valueEulDHShort    = eulerDH(N + 1) - strike - 0.5;
       valueEulDHAVShort  = eulerDHAV(N + 1) - strike - 0.5;
       valueEulDHSShort   = eulerDHS(N + 1) - strike - 0.5;
       valueEulDHLong     = eulerDH(N + 1) - strike + 0.5;
       valueEulDHAVLong   = eulerDHAV(N + 1) - strike + 0.5;
       valueEulDHSLong    = eulerDHS(N + 1) - strike + 0.5;
 
       % Compute the lower and upper bound for the nominator of the delta
       % hedge.
       payOffApproxInt(j)     = -max(valueIntShort, 0) + ...
         max(valueIntLong, 0);
         
       payOffApproxIntAV(j)   = 0.5 * (payOffApproxInt(j) - ...
         max(valueIntAVShort, 0) + max(valueIntAVLong, 0));

  	   payOffDHApproxInt(j)     = -max(valueIntDHShort, 0) + ...
         max(valueIntDHLong, 0);
         
       payOffDHApproxIntAV(j)   = 0.5 * (payOffDHApproxInt(j) - ...
         max(valueIntDHAVShort, 0) + max(valueIntDHAVLong, 0));
       
       payOffDHSApproxInt(j)  = -max(valueIntDHSShort, 0) + ...
         max(valueIntDHSLong, 0);

       % Compute all the necessary values for the Euler approximation.
 	   payOffApproxEul(j)     = -max(valueEulShort, 0) + ...
         max(valueEulLong, 0);
         
       payOffApproxEulAV(j)   = 0.5 * (payOffApproxEul(j) - ...
         max(valueEulAVShort, 0) + max(valueEulAVLong, 0));

  	   payOffDHApproxEul(j)     = -max(valueEulDHShort, 0) + ...
         max(valueEulDHLong, 0);
         
       payOffDHApproxEulAV(j)   = 0.5 * (payOffDHApproxEul(j) - ...
         max(valueEulDHAVShort, 0) + max(valueEulDHAVLong, 0));

       payOffDHSApproxEul(j)  = -max(valueEulDHSShort, 0) + ...
         max(valueEulDHSLong, 0);
    end
       
   % Compute the delta hedge
   deltaHedgeIntegral(iter)    = exp(-interestRate * T ) * ...
      sum(payOffDHApproxInt - payOffApproxInt) / (epsilon * M);
   deltaHedgeIntegralAV(iter)    = exp(-interestRate * T ) * ...
      sum(payOffDHApproxIntAV - payOffApproxIntAV) / (epsilon * M);
   % WARNING: The value bellow are not to be trusted!
   % \Delta = (V(S_0 + epsilon) - 2 * V(S_0) + V(S_0 - epsilon)) / 
   %   (2 * epsilon)
   deltaHedgeIntegralD(iter)   = exp(-interestRate * T) * ...
      sum(payOffDHApproxInt - 2 * payOffApproxInt + payOffDHSApproxInt) ...
      / (2 * epsilon * M);

   deltaHedgeEuler(iter)    = exp(-interestRate * T ) * ...
      sum(payOffDHApproxEul - payOffApproxEul) / (epsilon * M);
   deltaHedgeEulerAV(iter)    = exp(-interestRate * T ) * ...
      sum(payOffDHApproxEulAV - payOffApproxEulAV) / (epsilon * M);
   % WARNING: The values bellow are not to be trusted!
   % \Delta = (V(S_0 + epsilon) - 2 * V(S_0) + V(S_0 - epsilon)) / 
   %   (2 * epsilon)
   deltaHedgeEulerD(iter)   = exp(-interestRate * T) * ...
      sum(payOffDHApproxEul - 2 * payOffApproxEul + payOffDHSApproxEul) ...
      / (2 * epsilon * M);

end

% vim: expandtab :ts=4:tw=80

